%% %% EMAS IRFs
% This has been tested in Dynare 5.1.

clear all
tic
%addpath('Z:\dynare_programs\dynare-4.6.3\matlab');

dynare bbeffectivedemandmatchirf_order3.mod

toc

y_pos       = strmatch('y',M_.endo_names,'exact');
c_pos       = strmatch('c',M_.endo_names,'exact');
inv_pos     = strmatch('inv',M_.endo_names,'exact');
mu_pos      = strmatch('mu',M_.endo_names,'exact');
n_pos       = strmatch('n',M_.endo_names,'exact');
w_pos       = strmatch('w',M_.endo_names,'exact');
pi_pos      = strmatch('pie',M_.endo_names,'exact');
r_pos       = strmatch('r',M_.endo_names,'exact');
R_R_pos     = strmatch('rr',M_.endo_names,'exact');
sigma_a_pos = strmatch('vola',M_.endo_names,'exact');
Z_pos       = strmatch('z',M_.endo_names,'exact');
conditional_variance_R_E_pos = strmatch('varexpre',M_.endo_names,'exact');


IRF_periods      = 30;
burnin           = 2000; %periods for convergence
options_.pruning = 1;
options_.order   = 3;

shock_mat_with_zeros    = zeros(burnin+IRF_periods,M_.exo_nbr); %shocks set to 0 to simulate without uncertainty
IRF_no_shock_mat        = simult_amd(M_,options_,oo_.dr.ys,oo_.dr,shock_mat_with_zeros,options_.order)'; %simulate series
stochastic_steady_state = IRF_no_shock_mat(1+burnin,:); % stochastic_steady_state/EMAS is any of the final points after burnin

shock_mat                = zeros(burnin+IRF_periods,M_.exo_nbr);
shock_mat(1+burnin,strmatch('evola',M_.exo_names,'exact'))       = 1;
IRF_mat                  = simult_amd(M_,options_,oo_.dr.ys,oo_.dr,shock_mat,options_.order)';
IRF_mat_percent_from_SSS = (IRF_mat(1+burnin+1:1+burnin+IRF_periods,:)-IRF_no_shock_mat(1+burnin+1:1+burnin+IRF_periods,:)); %only valid for variables not yet logged

shock_mat(1+1+burnin,strmatch('evola',M_.exo_names,'exact'))     = 1;
IRF_mat                    = simult_amd(M_,options_,oo_.dr.ys,oo_.dr,shock_mat,options_.order)';
IRF_mat_percent_from_SSS_2 = (IRF_mat(1+burnin+1:1+burnin+IRF_periods,:)-IRF_no_shock_mat(1+burnin+1:1+burnin+IRF_periods,:)); %only valid for variables not yet logged


shock_mat(1+1+1+burnin,strmatch('evola',M_.exo_names,'exact'))   = 1;
IRF_mat                    = simult_amd(M_,options_,oo_.dr.ys,oo_.dr,shock_mat,options_.order)';
IRF_mat_percent_from_SSS_3 = (IRF_mat(1+burnin+1:1+burnin+IRF_periods,:)-IRF_no_shock_mat(1+burnin+1:1+burnin+IRF_periods,:)); %only valid for variables not yet logged

shock_mat(1+1+1+1+burnin,strmatch('evola',M_.exo_names,'exact')) = 1;
IRF_mat                    = simult_amd(M_,options_,oo_.dr.ys,oo_.dr,shock_mat,options_.order)';
IRF_mat_percent_from_SSS_4 = (IRF_mat(1+burnin+1:1+burnin+IRF_periods,:)-IRF_no_shock_mat(1+burnin+1:1+burnin+IRF_periods,:)); %only valid for variables not yet logged


irfmat_1_0 = IRF_mat_percent_from_SSS;
irfmat_2_1 = IRF_mat_percent_from_SSS_2 - IRF_mat_percent_from_SSS;
irfmat_3_2 = IRF_mat_percent_from_SSS_3 - IRF_mat_percent_from_SSS_2;
irfmat_4_3 = IRF_mat_percent_from_SSS_4 - IRF_mat_percent_from_SSS_3;


relv = strmatch('y',M_.endo_names,'exact');

figure(1);

plot(irfmat_1_0(1:end-3,relv),'Color',[0 0 0],'Linewidth',2);                           hold on;
plot(irfmat_2_1(2:end-2,relv),'Color',[0.2 0.2 0.2],'Linewidth',2,'Linestyle','--');    hold on;
plot(irfmat_3_2(3:end-1,relv),'Color',[0.4 0.4 0.4],'Linewidth',2,'Linestyle',':');     hold on;
plot(irfmat_4_3(4:end,relv),  'Color',[0.6 0.6 0.6],'Linewidth',2,'Linestyle','-.');    hold on;
plot(0*irfmat_4_3(:,relv),    'Color',[0 0 0]);                                         hold on;
title(['Order = ' num2str(options_.order) ', Pruning = ' num2str(options_.pruning)]);
xlim([0 20]);
ylim([-0.004 0.0005]);

stop;  %% Run code below to generate GIRF instead of EMAS IRF

%% %% GIRF
options_.pruning = 1;
options_.order   = 3;
%%eliminate shocks with 0 variance
i_exo_var        = setdiff([1:M_.exo_nbr],find(diag(M_.Sigma_e) == 0 )); %finds shocks with 0 variance
nxs              = length(i_exo_var); %number of those shocks
chol_S           = chol(M_.Sigma_e(i_exo_var,i_exo_var));%get Cholesky of covariance matrix to generate random numbers

shock_mat_with_zeros    = zeros(burnin+IRF_periods,M_.exo_nbr); %shocks set to 0 to simulate without uncertainty
IRF_no_shock_mat        = simult_amd(M_,options_,oo_.dr.ys,oo_.dr,shock_mat_with_zeros,options_.order)'; %simulate series
stochastic_steady_state = IRF_no_shock_mat(1+burnin,:); % stochastic_steady_state/EMAS is any of the final points after burnin
starting_point          = stochastic_steady_state; %define starting point of simulations; here: start at steady state and use 200 periods burnin to get to ergodic distribution

impulse_vec             = zeros(1,M_.exo_nbr); %initialize impulse vector to 0
impulse_vec(strmatch('evola',M_.exo_names,'exact')) = 1; %impulse of 1 percent of GDP shock to g

tic
shocks_baseline = [];
drop_periods    = 2000;
irf_replication = 1000; %Set to 2500 or lower number so that it computes faster

    IRF_mat = []; IRF_mat2 = []; IRF_mat3 = []; IRF_mat4 = [];
    irf_periods = 20;
    IRF_mat     = NaN(M_.endo_nbr,irf_periods,irf_replication);  %initialize matrix
    IRF_mat2    = NaN(M_.endo_nbr,irf_periods,irf_replication); %initialize matrix
    IRF_mat3    = NaN(M_.endo_nbr,irf_periods,irf_replication); %initialize matrix
    IRF_mat4    = NaN(M_.endo_nbr,irf_periods,irf_replication); %initialize matrix

for irf_iter = 1: irf_replication
    irf_iter
    shocks_baseline(:,i_exo_var)      = randn(irf_periods+drop_periods,nxs)*chol_S; %generate baseline shocks
    shocks_impulse                    = shocks_baseline; %use same shocks in impulse simulation
    shocks_impulse(drop_periods+1,:)  = shocks_impulse(drop_periods+1,:)+impulse_vec; %add deterministic impulse
    shocks_impulse2                   = shocks_impulse;
    shocks_impulse2(drop_periods+2,:) = shocks_impulse2(drop_periods+2,:)+impulse_vec; %add deterministic impulse
    shocks_impulse3                   = shocks_impulse2;
    shocks_impulse3(drop_periods+3,:) = shocks_impulse3(drop_periods+3,:)+impulse_vec; %add deterministic impulse
    shocks_impulse4                   = shocks_impulse3;
    shocks_impulse4(drop_periods+4,:) = shocks_impulse4(drop_periods+4,:)+impulse_vec; %add deterministic impulse
    
    y_baseline=[]; y_shock=[]; y_shock2=[]; y_shock3=[]; y_shock4=[];
    y_baseline = simult_amd(M_,options_,starting_point',oo_.dr,shocks_baseline,options_.order); %baseline simulation
    y_shock    = simult_amd(M_,options_,starting_point',oo_.dr,shocks_impulse,options_.order); %simulation with shock
    y_shock2   = simult_amd(M_,options_,starting_point',oo_.dr,shocks_impulse2,options_.order); %simulation with shock
    y_shock3   = simult_amd(M_,options_,starting_point',oo_.dr,shocks_impulse3,options_.order); %simulation with shock
    y_shock4   = simult_amd(M_,options_,starting_point',oo_.dr,shocks_impulse4,options_.order); %simulation with shock

    IRF_mat(:,:,irf_iter)  = (y_shock(:,M_.maximum_lag+drop_periods+1:end)-y_baseline(:,M_.maximum_lag+drop_periods+1:end)); %add up difference between series
    IRF_mat2(:,:,irf_iter) = (y_shock2(:,M_.maximum_lag+drop_periods+1:end)-y_shock(:,M_.maximum_lag+drop_periods+1:end)); %add up difference between series
    IRF_mat3(:,:,irf_iter) = (y_shock3(:,M_.maximum_lag+drop_periods+1:end)-y_shock2(:,M_.maximum_lag+drop_periods+1:end)); %add up difference between series
    IRF_mat4(:,:,irf_iter) = (y_shock4(:,M_.maximum_lag+drop_periods+1:end)-y_shock3(:,M_.maximum_lag+drop_periods+1:end)); %add up difference between series
    
    tim=toc;
    tim/irf_iter
    tim/irf_iter*(irf_replication-irf_iter)/60
end

GIRF_positive=[]; GIRF_positive2=[]; GIRF_positive3=[]; GIRF_positive4=[];
GIRF_positive  = mean(IRF_mat,3); %take average
GIRF_positive2 = mean(IRF_mat2,3); %take average
GIRF_positive3 = mean(IRF_mat3,3); %take average
GIRF_positive4 = mean(IRF_mat4,3); %take average

relv=26;
figure(2);
plot([0 GIRF_positive(relv,1:end-3)],'Color',[0 0 0],'Linewidth',2); hold on;
plot([0 GIRF_positive2(relv,2:end-2)],'Color',[0.2 0.2 0.2],'Linewidth',2,'Linestyle','--'); hold on;
plot([0 GIRF_positive3(relv,3:end-1)],'Color',[0.4 0.4 0.4],'Linewidth',2,'Linestyle',':'); hold on;
plot([0 GIRF_positive4(relv,4:end)],'Color',[0.6 0.6 0.6],'Linewidth',2,'Linestyle','-.'); hold on;
plot(0*GIRF_positive(relv,:),'Color',[0 0 0]); hold on;
title(['Order = ' num2str(options_.order) ', Pruning = ' num2str(options_.pruning)]);
xlim([0 20]);
ylim([-0.004 0.0005]);